# include "scene.h"

#ifdef SOLVE_ENERGY

# include "discr.h"
# include <stdio.h>

#ifdef _OMP_TEMP_ADVANCE_
# include <omp.h>
#endif

/*
# include <time.h>
# include <sys/timeb.h>
# include <stdlib.h>*/

# ifdef THMCOND_TDEP
extern double thm_cond(double T);
# endif

extern void temp_bc(const struct grid *g, double ***T, short tflag, unsigned short iRK/*, double *pE_accum*/);

/* Solving the Energy Equation using 4-stage RK-3 scheme. */
extern double temp_advance(struct grid *g)
  {
   short i, j, k, nr, nth, nz;
   double ***_T[5], ***u, ***v, ***w, ***T, ***T0, *r, *th, *z, density0, E_absr;
   char vhsrc_flag;
#  ifndef  THMCOND_TDEP
   double thmcond0;
#  endif

// struct timeb t1, t2;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'temp_advance()...'\n");
#  endif

   density0 = g->density0;
#  ifndef  THMCOND_TDEP
   thmcond0 = g->thmcond0;
#  endif

   _T[0] = g->T;
   _T[1] = g->arr.q[7]; _T[3] = g->arr.q[7];
   _T[2] = g->arr.q[8]; _T[4] = g->arr.q[8];

   T0 = g->T;
   u = g->um; v = g->vm; w = g->wm; r = g->r; th = g->th; z = g->z;
   nr = g->nr; nth = g->nth; nz = g->nz;


   if(g->vhsrc_flag == 'y') if(T0[NrG+nr/3][NthG+nth/2][NzG]-T_INLET > EXPECTED_T_RISE) g->vhsrc_flag = 'n';
   vhsrc_flag = g->vhsrc_flag;

// ftime(&t1);
// printf("-------\n");

/* ------------------- Calculating T(1) ------------------- */
   T = _T[0];
#  ifdef _OMP_TEMP_ADVANCE_
   #pragma omp parallel for private(i,j,k) schedule(guided)
#  endif
   for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
   for(k = NzG; k < nz + NzG; k++) _T[1][i][j][k] = _T[0][i][j][k] + 0.5*g->dt*E_ijk;
// ftime(&t2);
// printf("dt = %ld ms\n", t2.time*1000+t2.millitm-t1.time*1000-t1.millitm);
// exit(-1);

   temp_bc(g, _T[1], 0, 1);
/* -------------------------------------------------------- */

/* ------------------- Calculating T(2) ------------------- */
   T = _T[1];
#  ifdef _OMP_TEMP_ADVANCE_
   #pragma omp parallel for private(i,j,k) schedule(guided)
#  endif
   for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
   for(k = NzG; k < nz + NzG; k++) _T[2][i][j][k] = _T[1][i][j][k] + 0.5*g->dt*E_ijk;

   temp_bc(g, _T[2], 0, 2/*, &pE_accum*/);
/* -------------------------------------------------------- */

/* ------------------- Calculating T(3) ------------------- */
   T = _T[2];
#  ifdef _OMP_TEMP_ADVANCE_
   #pragma omp parallel for private(i,j,k) schedule(guided)
#  endif
   for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
   for(k = NzG; k < nz + NzG; k++) _T[3][i][j][k] = (2.0/3.0)*_T[0][i][j][k] + (1.0/3.0)*_T[2][i][j][k] + (1.0/6.0)*g->dt*E_ijk;

   temp_bc(g, _T[3], 0, 3/*, &pE_accum*/);
/* -------------------------------------------------------- */

/* ------- Getting "T" at next time-step. It will be given by "_T[4]" --------- */
   T = _T[3];
#  ifdef _OMP_TEMP_ADVANCE_
   #pragma omp parallel for private(i,j,k) schedule(guided)
#  endif
   for(i = NrG; i < nr + NrG; i++) for(j = NthG; j < nth + NthG; j++)
   for(k = NzG; k < nz + NzG; k++) _T[4][i][j][k] = _T[3][i][j][k] + 0.5*g->dt*E_ijk;

   temp_bc(g, _T[4], 1, 0/*, &pE_accum*/);
/* ------------------------------------------------------------------------------ */


/* Determining total energy absorbed by the system. */
   E_absr = 0;
#  ifdef _OMP_TEMP_ADVANCE_
   #pragma omp parallel private(i,j,k)
#  endif
     {
#     ifdef _OMP_TEMP_ADVANCE_
      #pragma omp for schedule(guided) reduction(+:E_absr) nowait
#     endif
      for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++) for(k = NzG+nz-1; k >= NzG; k--)
        E_absr += (r[i+1]+r[i])*(r[i+1]-r[i])*(th[j+1]-th[j])*(z[k+1]-z[k])*(_T[4][i][j][k] - _T[0][i][j][k]);

#     ifdef _OMP_TEMP_ADVANCE_
      #pragma omp for schedule(guided)
#     endif
      for(i = 0; i < nr + 2*NrG; i++) for(j = 0; j < nth + 2*NthG; j++) for(k = 0; k < nz + 2*NzG; k++)
   /* Copying n+1-th step temperature to _T[0]. */
       _T[0][i][j][k] = _T[4][i][j][k];
     }

   E_absr = E_absr*0.5*g->density0*Cp_WATER;

   g->T = _T[0];

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'temp_advance()' ends.\n");
#  endif

   return E_absr;
  }
#endif
